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Abstract 

We have examined, repeated and extended earlier numerical calculations of Berger and 
Moncrief for the evolution of unpolarized Gowdy T 3 cosmological models. Our results are 
consistent with theirs and we support their claim that the models exhibit AVTD behaviour, 
even though spatial derivatives cannot be neglected. The behaviour of the curvature invariants 
and the formation of structure through evolution both backwards and forwards in time is 
discussed. 
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1 Introduction 

The Hawking-Penrose theorems || imply that every 'physically reasonable' cosmological model 
possesses a singularity; however very little is known about its nature. Early studies due to Belinskii 
et al [0 suggested that the behaviour of the generic singularity resembles that of the spatially 
homogeneous Bianchi VIII and IX cosmological models, commonly called mixmaster dynamics. 
It is generally felt, however, that matters are not this simple, and attention has focused on the 
velocity-term-dominated behaviour (VTDB) hypothesis due to Eardley et al jjj. It is very hard to 
give a covariant definition but, loosely speaking, the hypothesis claims that near the singularity 
spatial derivatives play no part in the evolution; locally one has a spatially homogeneous cosmology, 
but with the parameters varying from point to point. This is the starting point for many large- 
scale-structure calculations in physical cosmology (see, for example, |lq]), but is it correct? If 
spatial derivatives are ignored, the evolution equations become ordinary differential equations 
whose behaviour is easier to analyse. Isenberg and Moncrief have introduced an alternative 
hypothesis of asymptotically velocity-term dominated (AVTD) behaviour. A cosmological model is 
AVTD if its asymptotic behaviour near the singularity is that of the ordinary differential equation 
system. Clearly VTDB implies AVTD behaviour, but the converse may not be true. In order to test 
the validity and usefulness of these hypotheses we need to examine the behaviour of inhomogeneous 
cosmological models, and the simplest examples would appear to be the Gowdy T 3 models ||. 
Polarized Gowdy cosmologies exhibit AVTD behaviour [JlO|, and so attention focuses on the less 
tractable unpolarized case M, which appears to require numerical treatment. 

Berger and Moncrief and Berger || (who give much more extensive bibliographies than 
the brief introduction above) have tackled the unpolarized case numerically. This is a decidedly 
non-trivial task, and so they introduced novel numerical algorithms, which led them to conclude, 
tentatively, that the AVTD hypothesis was confirmed. Their grounds for caution are based on 
the inability of their code to resolve fine-scale spatial structure. As they suggest, this problem 
requires an algorithm with adaptive mesh refinement (AMR), discussed in section 3. 



We have used our AMR code to repeat and extend the earlier calculations. We obtain broadly 
similar results, but with additional fine structure. We confirm the earlier result that the Gowdy T 3 
cosmologies exhibit AVTD behaviour but can find no evidence that the VTDB hypothesis holds 
universally near the singularity. Our results show that by starting with smooth data and evolving 
backwards in cosmic time towards the singularity our fields develop complicated spatial structure. 
If we time-reverse our calculations we appear to be asserting that complicated structure becomes 
smooth as it evolves into the future, the exact opposite of widespread beliefs in the theory of large- 
scale structure! We have therefore taken generic smooth data near the singularity and evolved it 
into the future finding, as expected, that complicated spatial structure appears. The fields that 
we are discussing are merely metric components with no covariant meaning. We have therefore 
computed curvature invariants. These become very large near the singularity, but they too have 
complicated spatial structure, implying that this is not a coordinate effect. 



2 The Gowdy T 3 Universe 

The line element used by Berger and Moncrief |J for the Gowdy T 3 cosmology is 
ds 2 = e A /V/ 2 (- e - 2r dr 2 + dd 2 ) + 

(2.1) 

e- T (e p da 2 + 2e p Q da dS + (e p Q 2 + e - p ) dd 2 ), 

where A, P and Q are functions of r and 9 only, — oo < t < oo with t = oo a singularity, and 
< 8, a, S < 2ir. The functions A, P and Q are required to be periodic in 8 with period 2ir. The 
polarized mode corresponds to Q = 0. 

The vacuum momentum constraint equation is 

-X e = 2(P T Pg + e 2P Q T Q e ), (2.2) 

and the Hamiltonian constraint is 

-K = Pr 2 + e- 2T Pe 2 + e 2P (Q T 2 + e^Qg 2 ), (2.3) 
while the vacuum evolution equations reduce to 

P TT = e- 2T P ee + e 2P {Q 2 - e^Qg 2 ), (2.4a) 

Qtt = e- 2r Qgg - 2(P T Q T - e- 2r PgQg). (2.4b) 

Here f T = df/dr, etc. It should be noted that the evolution equations ( |2.4|) do not involve A. 



The energy constraint (2.3) determines A and the momentum constraint (2.2) is satisfied at all 
times if it is satisfied initially, provided the other equations hold. As pointed out by Berger and 
Moncrief Q| the evolution equations are harmonic map equations for a target space with metric 

dS 2 = dP 2 + e 2P dQ 2 , (2.5) 

and so study of these equations is not without interest. 

The integration programmes described in the next section require a first order system of equa- 
tions. We accomplish this by introducing new variables as follows: 

A = P T , B = Q T , C = P e , D = Q e . (2.6) 

The system of evolution equations now takes the form 

(2.7a) 
(2.7b) 
(2.7c) 
(2.7d) 
(2.7e) 
(2.7f) 

2T D 2 ). (2.7g) 
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The last equation is only needed if A is to be evolved simultaneously with P and Q, in which case 

\ e = -2{AC + e 2P BD) (2.8) 

from equation ( |2 . 2[ ) acts as a constraint. 

Initial data at r = is also required. Berger and Moncrief suggest that the following choice is 
reasonably generic (with A chosen to satisfy fl2.8p ) : 

A = wo cos 0, 5 = 0, (7 = 0, £> = -sin0, P = 0, Q = cos0, A = 0, (2.9) 

and the results reported here use vq = 10 as used in ||. 



3 Numerical Procedures 

As we shall see in the next section, the evolution produces structure on fine scales with steep 
gradients. One numerical strategy is to use a standard algorithm, e.g., Lax-Wendroff, together 
with a very fine grid structure. This certainly works but it is not the most efficient way to proceed. 
Berger and Moncrief || used a symplectic integrator, described in their paper. However their 
calculation produced fine-scale structure even at sizes comparable with the grid spacing, which 
they regarded as unreliable. In presenting their results spatial averaging was used to remove the 
finest-scale structure. An alternative approach is due to van Putten Jl4|] , but he reports integration 
of the equations only for short times and small values of the parameter vq in the initial data. 

Problems of this type are natural candidates for adaptive mesh refinement (AMR), which seeks 
to add extra grid points when they are needed and to remove them when they become unnecessary. 
We used an implementation of the Berger and Oliger Q AMR algorithm which is descended from 
the one described in Hamade and Stewart ||. Numerous changes have been made to that code. 
Those relevant to this paper are: 

• an arbitrary number of grids may exist at each level of refinement; 

• it is now possible to switch in arbitrary integrators, giving the code 'modularity'; 

• cubic spline interpolation is preferred to the quadratic interpolation used earlier, thus guar- 
anteeing continuous spatial second derivatives; 

• when creating finer grids one can choose to use data from obsolete grids at the same level, 
rather than relying on interpolation from the parent grid. 

For non-specialists, the third change introduces a subtle form of data smoothing, which can be 
(almost) eliminated by the fourth change. 

It is straightforward to write a generalization of the standard Lax-Wendroff two-step algorithm 



which is second-order accurate for the system (2.7), and this proved to be the fastest choice. (If 



the typical spatial grid spacing is A9, then a second-order accurate method produces a local 
truncation error which is 0((A8) 3 ). The mesh refinement is typically AO —> A0/4, and so a single 
stage of refinement reduces the local truncation error by nearly two orders of magnitude, as well 
as enhancing spatial resolution.) In addition (and as part of a long-term programme) we used 
our interpretation of the second-order accurate wave propagation routine from the CLAWPACK 
high-resolution package of LeVeque [Q. (CLAWPACK is written in FORTRAN. The administrative 
part of our code is written in CH — h, and the numerically intensive functions are written either 
in C or C++. It proved easier, and more instructive, to adapt the CLAWPACK routines to this 
environment, rather than to interface the FORTRAN to the existing code.) 

For calculations using the Berger-Oliger algorithm the initial coarse grid had a spatial sep- 
aration A0 = 2?r/2000. The finest child grid actually used had AO = 2tt/512 000. As a check 
to ensure that the AMR code was not introducing spurious effects we repeated the calculations 
with Ad = 27t/8000 and no mesh refinement. The Lax-Wendroff and CLAWPACK- like codes were 



3 



T = QlT 



T = 547T 





-500 



0.5 



1.5 



2.5 



-500 1 



0.5 



2.5 




Figure 1: P, Q and A as functions of 6 at two different times. The left column corresponds 
to t — 6-7T, the right to r = 54-7T, and 9 ranges from to it. (The problem is invariant under 
9 — > 2-7T — 8, and so only one half of the full 6-range is shown.) The large spikes in Q are 
truncated by the choice of vertical scale. 



developed by the two different authors, continuing our policy to eliminate coding errors. We also 
checked our codes against the exact pseudo-unpolarized solution 

P = log cosh P , Q = tanhP , P (t, 9) = J {e' T ) cos 6, (3.1) 

(where Jo is a Bessel function) given by Berger and Moncrief , and confirmed that for both of our 
integration methods the errors decreased at second order. There were only very small differences 
between the results from AMR with Lax-Wendroff and those from AMR with CLAWPACK, and 
what is reported in the next section is common to both. 



4 Numerical Results 

We are reporting on the integration of equations Jp.Tj ) with initial data ( |2.9D and vq = 10, the 
same conditions as used by Berger and Moncrief || . (More or less spatial structure develops for 
larger or smaller, respectively, values of the parameter v . Berger g] uses v = 5 and gets results 
with a less complicated spatial structure, though with the same basic features.) 

The main numerical results of Berger and Moncrief are summarized in figures 4 and 5 of |J . 
The left-hand column of our figure 1 shows P, Q and A as functions of 9 at r = 6ir = 18.85, the 
latest time of their figures 4 and 5. The right-hand column of figure 1 shows the variables at a 
much later time of r = 547T = 169.6. The graphs in figure 1 comprise all available data from the 
AMR grid hierarchy. 

The structure of P becomes complicated soon after the start of the evolution and by time 
t = 67r (as figure 1 shows) it includes a number of fine spikes. At first sight these spikes might 
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Figure 2: Data from three generations of grids at a spike in P. Each grid has a resolution 
greater by a factor of four than the grid to its left. Each circle is one grid point, and the 
dotted boxes show the sizes and data ranges of the subgrids. (The data is taken at time 
r = 15.) 



seem to be errors in the evolution, and, indeed, when Berger and Moncrief first found spikes 
in their data, they considered them unreliable and used spatial averaging to remove them from 
figures 4 and 5 of 0. However, our AMR code automatically refines the evolution in regions 
where spikes start to form, increasing the grid resolution so that the spikes appear smooth and 
the evolution can continue without loss of accuracy. This is demonstrated in figure 2 where data 
at a spike is shown for a coarse grid and for two refined subgrids. (Notice that on the finer grids 
the spike reaches a greater height than on the coarse grid. The AMR algorithm uses an average 
of data from fine grids to update data on coarse grids.) 

Spikes form in P in both the positive and negative directions. Examining the development 
of the positive spikes (the negative spikes are discussed below) it is found that their peaks grow 
linearly at rates P T = constant > 1, while at the same time they decrease in width. As a spike 
narrows (and there is no indication that the spikes ever stop growing or narrowing) the AMR code 
is forced to use finer and finer grids to keep it well resolved until, eventually, it becomes narrower 
than the smallest grid spacing that the AMR code has been instructed to use, and from that point 
on it cannot be assumed that the code is accurately evolving the spike. We have to accept that 
there are probably significant errors in the calculated heights of the spikes in P, but this does not 
imply that there must be large errors elsewhere in the data: since the characteristic speed of the 
model is exp(— r), deviations from the exact solution that develop at some time not too soon after 
the start of the evolution can only propagate very small distances from their points of origin. 

At late times P TT = A T becomes very small, and so A freezes and P grows linearly in time 
with the spatial profile shown in the right-hand column of figure 1. Except at the points where 
spikes formed, the spatial profile of P at late times has quite a simple form, much less intricate 
than early on in its evolution. For reasons discussed above, the numerical simulation is unable 
to determine the late-time behaviour of the spikes. While small spikes and lumps appear in the 
graph of P at t = 547r in figure 1, these are just the residue left by spikes that have narrowed 
beyond the resolution of the simulation. The spikes themselves, assuming they continued to grow 
at their initial rates, would go off the scale of the graph. 

To suggest reasons for the observed behaviour it is useful to know which are the dominant terms 
in the governing equations. Here we adopt the definition that a term is dominant if its absolute 
value exceeds ten times the sum of the absolute values of the other terms. Figure 3 illustrates 
this for the wave equations (2.7a,b) for A T = P TT a nd B T = Q TT . Th e lig htly shaded regions in 
(r, (9)-space are where the terms exp(2P)f? 2 (in (2.7a)) and —2AB (in (2.7b)) dominate the right- 
hand sides. The darker shaded regions are where the terms — exp(2P — 2t)D 2 and 2 exp(— 2t)CD 
dominate. The narrow black regions are where the terms exp(— 2r)Cg and exp(— 2r)Dg dominate. 
The white regions are where no one term dominates. It is clear from the figure that there is very 
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Figure 3: The dominant terms in the A T and B T equations (2.7 a,b) as functions of 6 and r. 
The black regions are where terms involving Cg and Dg dominate, the dark-shaded regions 
are where terms involving C and D dominate, and the light-shaded regions are where terms 
involving A and B dominate. In the white regions no single term dominates. 



complicated interplay between terms during the early stages of the evolution, and it seems unlikely 
that the overall behaviour can be traced to simple causes. Berger || notes that spikes form at 
points where D = Qg = (though not all such zeros produce spikes) and discusses their formation 
in terms of bounces off of potentials in the evolution equations. In terms of dominance behaviour, 
spike formation occurs in dark-grey regions in the left half, and light-grey regions in the right half 
of figure 3. 

Sharp features also form in Q as it evolves. Figure 1 shows that Q at r = 6tt has a negative 
spike and a positive spike very close together with a steep gradient between the two peaks. A 
second, similar double-spike feature is also present, but only just visible in figure 1 because of the 
scale. The positions of the two double spikes in Q coincide with the positions of the two negative 
spikes in P. A close examination of the data reveals that both the positive and negative peaks of 
each double spike in Q grow at an exponential rate, while the overall width of the feature becomes 
smaller, and that the negative spikes in P grow linearly in much the same way that their positive 
counterparts do. The discussion above regarding the numerical problems caused by the narrowing 
of spikes applies equally here. 

A double-spike feature is also reported in ^ (see especially figure 2 there, but note that it is 
produced using different initial data) , however the explanation given there is not quite convincing. 
Examination of the dominance data of figure 3, together with similar data for finer grids, reveals 
that the region around a double spike (6 « 2, r > 4 for the large double spike in figure 1) is light 
grey for both equations. We may conclude that in this region the terms proportional to exp(— 2t) 
are irrelevant and so 



e 2P B\ 



B T ~ -2AB. 



(4.1) 



A double spike then forms under the conditions of B crossing zero (the exact position of this 
shows up as a thin white line in both halves of figure 3), while A is negative and P is negative or 
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small and positive. It follows that A T > and is small, so that A remains negative and P forms a 
negative spike. On both sides of the zero, B then grows exponentially, but in opposite directions, 
and the result is a double spike in Q. Since A T > the region in which A < becomes smaller, 
and the double spike narrows. (By r = 10 the light-grey region that supports the growth of the 
larger double spike is no longer visible in the left half of figure 3.) 

Away from the double spikes, Q T = B quickly approaches zero and Q freezes not long after 
the start of the evolution: the two plots of Q in figure 1 differ only at the double spikes. Because 
of the narrowing effect, the second double spike is still barely visible on the right half of figure 1 
although, extrapolating its initial growth, it should be larger than the vertical scale of the graph. 

The evolution of A shows similar features to the evolution of P (see figure 1). Negative spikes 
develop in A at the same points at which positive spikes develop in P and, apart from pointing in 
a different direction, they show the same behaviour: they grow at a constant rate and they narrow 
until they no longer can be resolved by the simulation. In contrast, no unusual behaviour seems 
to develop in A at the points where double spikes form in Q (although the errors in the data at 
the double spikes are sufficient to produce erratic behaviour in A there at late times, as can be 
seen in figure 1). 

Equation (EJg) shows that A is a decreasing function of time, and it is clear from figure 1 



that A becomes large and negative as it evolves. At late times, A T fss —A 2 (cf. equation ( 2.7 g)) 
and A is approximately constant, so A grows linearly with a fixed spatial profile. 



The constraint equation (2.8) in our system relates the spatial derivative of A to the values 



of P and Q and their first derivatives. For an exact solution of the evolution equations (2/7) the 
constraint equation is satisfied at all times if it is satisfied initially, but for a numerical solution 
this is not necessarily the case. We monitor the degree to which the constraint equation fails to 
hold by calculating the total percentage error defined by 

error(%) = 100 ^ \\ e + 2{AC + e 2P B D)\ d6 j J \\ e \d6, (4.2) 

where the integrals are evaluated using the trapezium rule and Xg is estimated from A using 
second-order finite differencing. Since the errors are known to be disproportionately large at the 
spikes we do not include these points when measuring the constraint. (A region of width 0.02 is 
excised at each spike, with about ten percent of the total domain being removed.) At time r = 6tt 
the error in the constraint is 1.1%. At time r = 54ir the error is 0.97%, smaller because of the 
overall growth of A. 



Instead of evolving A with equation ( 2.7 g) we could calculate it at every time step by integrating 
equation (2.8). This would obviously ensure that the constraint equation is satisfied, and such 
'fully constrained' approaches are often applied to problems in numerical relativity. However, for 
this particular problem enforcing the constraint equation in this way would be a very bad idea: 
the large errors at the spikes would affect the whole of A rather than being confined to isolated 
regions. 

Our understanding of a spacetime is, in general, advanced more by studying coordinate- 
independent quantities which have some physical relevance than by studying the components of 
the metric. For a vacuum spacetime we can use the curvature tensor to construct four independent 
scalar quantities — the curvature invariants — which we take in the form 



,„ T _ }_n fifivpo _n p<?n *? w f7 ^ v 

Wj — — O^j/po-O , win — — o p(J , 



(4.3) 



where C pvprT is the Weyl tensor, e pvpa is the Levi-Civita tensor and C* vpa = \s ^^C^ p<y . 

Figure 4 shows two of the curvature invariants, w\ and wu, evaluated using our data at time 
t = 6tt. The following discussion concentrates on these two quantities since wm and wiy behave 
qualitatively the same as, respectively, wi and wu. 
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Figure 4: Curvature invariants W\ and W\\ (see equations (4.3)) at r = Qtt for < 9 < tt 



While w\ is invariant under the reflection 9 — > 2ty — 9, the sign of tun is reversed. (The data 
has been scaled using the transformation w — ► (1/ In 10) sinh _1 («;/2).) 



The most striking feature of the curvature invariants is their magnitude: Wj and u>n at r = 6n 
reach both positive and negative values of size greater than ten to the hundredth power. At the 
start of the evolution w\ takes values from —50 to +30 000, and wu from —2500 to +2500. Their 
subsequent growth is rapid, and shows no signs of reversing. The main influence on the size of the 
curvature invariants is the factor exp(— A) in the expressions for wi and wu. The effect of A can 
be seen by comparing the profiles of W\ in figure 4 and A at r = 67r in figure 1. The linear growth 
(towards minus infinity) of A at late times corresponds to an exponential growth of the curvature 
invariants. 

As figure 4 shows, the curvature invariants change sign suddenly at various values of 9. These 
sign changes tend to form into clusters which contract in width as the evolution progresses. In w\ 
the clusters consist of even numbers of sign changes and the function is predominantly positive, 
while in w\\ they are odd numbered and the function changes sign between clusters. The sign 
change clusters become narrow in just the same way that the spikes in P do, and indeed the 
two types of feature coincide. The narrowing, together with the inaccuracies in the data at the 
spikes, makes determining the precise behaviour of sign changes within a cluster difficult. However, 
careful examination suggests that while the width of a cluster tends to zero, there is no cancellation 
between the sign changes. 

Given the large numbers of terms in the expressions for the curvature invariants, relating their 
behaviour to features of P, Q and A is not a simple matter. However, by comparing the magnitudes 
of individual terms we have found that at late times we can approximate the curvature invariants 
with simplified expressions 

™i«^e 3T -V 2 +3)(^ 2 -l) 2 , (4.4a) 

wu « -—e 2T+p - x DA(A + 3)(A- l)(A + if. (4.4b) 
16 

These approximations fail when the values they predict are small: other terms in the curvature 
invariants clearly cannot be neglected then. The sign changes in w\ are all clustered around the 
positive spikes in P, and this makes sense in terms of the approximation since — 1<A<+Iis 
found everywhere except at these spikes. (The terms neglected in (4.4a) turn out to be negative 
at the spikes.) The behaviour of w\\ is similar, except for the influence of D in equation (4.4b). At 
each positive spike in P there is a zero of D together with an odd number of sign changes in wn. 
Zeros of D at points where no spikes develop produce isolated sign changes in However, 
although a pair of zeros appears in D where a double spike forms in Q, no unusual behaviour 
occurs in wu there. It seems that zeros in A counteract the effect of the zeros in D sufficiently so 
that terms neglected in equation (4.4b) become significant and prevent wu from changing sign. 

(We note that the sign change clusters in wj and wu evolve at different rates, and this leads 
to some discrepancies between figure 4 and the behaviour described in the paragraph above. By 
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Figure 5: Developmen t of small-scale structure in P for evolution both forwards and backwards 
in time. Initial data (2.9) with Vq — 15 is used at starting time r = 1.5. The left two plots 



are produced by evolving backwards from this initial state, the right two plots by evolving 
forwards. 



time t — 6ir some of the clusters are already too narrow to be resolved and inaccuracies at the 
double spikes are large enough to cause the invariants to behave erratically there.) 

It follows that the formation of positive spikes in P is more than a coordinate effect induced 
by the form of the Gowdy metric, and is connected with unusual behaviour of the spacetime 
curvature. On the other hand our results suggest that this may not be the case for the double 
spikes that form in Q: while both P and Q can be related to the proper lengths of the trajectories 
of the metric's T 2 isometry group, the spiky behaviour of Q does not seem to have a physical 
effect beyond this. 

If we extrapolate the observed behaviour of the curvature invariants forward in time to the 
singularity at r = oo we see that they become infinite — except possibly at isolated values of 6 
corresponding to the positive spikes in P. This is in agreement with the curvature behaviour 
at the singularity conjectured by Moncrief Jl3| . (The same reference shows that no curvature 
singularities can develop before r = oo, and our results do not suggest otherwise.) 

The main features of the behaviour reported above are also found in simulations using initial 
data other than equation (2.9) or a starting time other than r = 0, and the implication is that the 
behaviour is reasonably generic. It is interesting to ask, though, what behaviour the spacetime 
must to have prior to the starting time in order that it arrives at the initial state we have chosen — 
that is, what happens if we start with initial data ( |2.9| ) but reverse the direction of time? We find 
that, in general, the variables P and Q have a simple wave motion away from the singularity, as 



equations (2.4) suggest, but that for a starting time closer to the singularity (e.g., r > 4) or a 



large value for the parameter vq in (2J5) small-scale spatial structure again starts to form. (This is 
illustrated in figure 5.) However, in this case the small-scale structure does not grow indefinitely 
and eventually the wave motion takes over. Even when small-scale structure is developing, the 
curvature invariants are found to be rapidly approaching zero, and the spacetime 'flattens' as it 
expands. 



5 The AVTD and VTDB Hypotheses 

Although the VTDB hypothesis of Eardley et al [|| is difficult to formul ate in a co vari ant manner, 
it is easy to interpret it for this problem. In our evolution equations (2J3) and (2^4) (combined 
in (I]?])) every spatial derivative is multiplied by a factor exp(— r). Since the singularity is given 
by r — > +oo, it might seem reasonable to ignore these terms when studying the singularity. The 
resulting system of ordinary differential equations can be solved analytically and one finds 



P(t,0) ~a o (0)T, Q(T,8)~q (8) 



(5.1) 



The numerical evidence in favour of this asymptotic behaviour is compelling: B = Q T and A T = 
P TT tend to zero as r increases. Thus we can confirm Berger and Moncrief 's assertion in [^| that 
the Gowdy T 3 cosmologies exhibit AVTD behaviour. 
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Figure 6: The VTDB region for r < 54vr. Spatial 
derivatives are considered negligible in the shaded 
regions, which correspond to the intersection of 
the light grey regions from both halves of figure 3. 
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However, there is no evidence that the VTDB hypothesis used to derive (5.1) is valid! The 
light-grey regions in figure 3 show where spatial derivative terms are unimportant in the two wave 
equations, and examination of their intersection reveals that for r < 6ir there are intervals of 
in which VTDB occurs, but that they do not cover the whole range of 9. Figure 6 shows these 
VTDB regions for an extended time range, r < 54-7T. Clearly the behaviour is not universal. 

In a recent preprint, Kichenassamy and Rendall [ [rT| have examined analytically the behaviour 
of the Gowdy T 3 co smologies near the singularity. Their results seem to suggest that the asymp- 
totic behaviour (5.1) is universal, and is therefore not contingent on the VTDB hypothesis. 



6 Conclusions 

We have two groups of conclusions, one for cosmology and one for numerical relativity. We have 
investigated structure formation and behaviour near the singularity in the Gowdy T 3 cosmologies. 

• Our calculations are consistent with the earlier ones of Berger and Moncrief g| , but show 
more fine-scale structure. 

• This fine-scale structure is not a coordinate effect since it also occurs in the curvature in- 
variants. 

• As claimed earlier there is considerable numerical evidence that the AVTD conjecture is 
valid for the Gowdy T 3 cosmologies, but there is no evidence to support the hypothesis of 
VTDB. Although spatial derivative terms matter, this does not seem to affect the asymptotic 
evolution near the singularity. 

• The Gowdy T 3 cosmologies have the unexpected property that smooth initial data acquires 
structure when it is evolved both backwards and forwards in time. 

We have not investigated in detail the symplectic integrator used by Berger and Moncrief, but 
for this problem it would seem to perform no better than standard finite-difference methods. A 
simple Lax-Wcndroff code using 8000 grid points without AMR will produce results that differ from 
those shown in figure 1 only at the points where spiky features develop. However, the AMR code 
is faster and allows details smaller than the coarse grid spacing to be investigated. Summarizing 
the advantages AMR brings for this problem: 
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• Results (on the coarsest grid) with the required level of overall accuracy are produced in 
minimal time. 

• Small spatial features are resolved without requiring a high density of grid points everywhere 
on the domain. 

• The total number of grid points changes with time and allows quick evolution of the 'frozen' 
spatial structures found at late times. 

• No prediction as to the amount of fine-scale structure that will develop is needed beforehand. 

There is another, subtle feature of AMR. We have been asked frequently why Lax-Wendroff 
with AMR gives results which are indistinguishable from the highly sophisticated CLAWPACK codes. 
At first sight the Lax-Wendroff scheme would appear to be inappropriate: it contains artificial 
viscosity which would smear out the spatial structure we observe. However the viscosity coefficient 
is 0((A9) 2 ) where A8 is the spatial grid spacing, and when spatial structure is encountered the 
AMR code evolves it with a very small A0 so that consequently there is very little dissipation. 

We thank Neil Cornish for useful discussions, and Beverly Berger, David Garfinkle and Vincent 
Moncrief for encouragement. Simon Hern was supported by an EPSRC studentship. 
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